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ABSTRACT 

The Cygnus region is a very bright and complex portion of the TeV sky, host to unidentified sources and a 
diffuse excess with respect to conventional cosmic-ray propagation models. Two of the brightest TeV sources, 
MGRO J2019+37 and MGRO J2031+41, are analyzed using Milagro data with a new technique, and their emission 
is tested under two different spectral assumptions: a power law and a power law with an exponential cutoff. 

The new analysis technique is based on an energy estimator that uses the fraction of photomultiplier tubes in 
the observatory that detect the extensive air shower. The photon spectrum is measured in the range 1-100 TeV 
using the last three years of Milagro data (2005-2008), with the detector in its final configuration. An /--test 
indicates that MGRO J2019+37 is better fit by a power law with an exponential cutoff than by a simple 
power law. The best-fitting parameters for the power law with exponential cutoff model are a normalization at 
10 TeV of T_\ x 10- 10 s” 1 m" 2 TeV" 1 , a spectral index of 2.0 40 ' 5 , and a cutoff energy of 29 + _f (j TeV. MGRO 
J2031+41 shows no evidence of a cutoff. The best-fitting parameters for a power law are a normalization of 2 .Dq 6 x 
10~ 10 s _1 m -2 TeV -1 and a spectral index of 3.224o 23 8 . The overall flux is subject to a ~30% systematic uncertainty. 

The systematic uncertainty on the power-law indices is ~0. 1 . Both uncertainties have been verified with cosmic-ray 
data. A comparison with previous results from TeV J2032+4130, MGRO J2031+41, and MGRO J2019+37 is also 
presented. 

Key words: acceleration of particles - astroparticle physics - gamma rays: general - open clusters and associations: 
individual (Cyg OBI, Cyg OB2) - pulsars: general 

Online-only material: color figures 


1. INTRODUCTION 

The Cygnus region is a part of our Galaxy of active massive 
star formation and destruction, and has been studied over a 
broad range of wavelengths, including radio, X-ray, GeV, and 
TeV gamma ray, as well as in cosmic rays. From GeV up to 
multi-TeV energies, the Cygnus region is the brightest diffuse 
gamma-ray source in the northern hemisphere (Hunter et al. 
1997). 
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One of the challenges in analyzing the Cygnus region at TeV 
energies is the proper separation of the gamma-ray flux that is 
attributed to the point or extended sources in the region or to 
the diffuse emission. Previous Milagro analyses computed the 
diffuse emission from the region using two alternative methods 
to isolate the contribution from the resolved sources (Abdo 
et al. 2007a, 2008), and found that at TeV energies the flux 
is still in excess with respect to the predicted flux from the 
GALPROP model (Strong & Moskalenko 1998; Moskalenko 
et al. 1998, 2000). Milagro also published the discovery of two 
unidentified TeV sources in the region, MGRO J20 19+37 and 
MGRO J2031+41 (Abdo et al. 2007b, 2009b). The location of 
MGRO J2019+37 was found to be consistent with two EGRET 
sources, while the best-fit position for MGRO J2031+41 was 
near two EGRET sources and the HEGRA unidentified source 
TeV J2032+4130. In a correlation study connecting the TeV 
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sources discovered by Milagro with sources detected above 10a 
(the so-called Bright Source List, BSL) by the Fermi Large Area 
Telescope (LAT), the two aforementioned brightest Milagro 
sources in the region were associated with GeV pulsars (Abdo 
et al. 2009a, 2009b). MGRO J20 19+37 and MGRO J2031+41 
are currently associated with two and one pulsars identified 
by the Fermi LAT, respectively (The Fermi LAT Collaboration 
2012b). Milagro also detected two candidate sources in the 
Cygnus region (/ e [65°, 85°] and b e [-2°, +2°]; Abdo et al. 
2009b). 

Recently, several collaborations have presented new surveys 
of the Cygnus region, confirming the complexity of the region 
and showing the highly structured diffuse emission. The TeV 
emission from the position of MGRO J2019+37 has been 
confirmed by the VERITAS experiment (Aliu et al. 2011). The 
VERITAS spatial counts map shows a clear structure associated 
with MGRO J2019+37, a region of extended emission which 
seems to be produced by previously unresolved sources. At 
lower energies, the Fermi LAT Collaboration (Ackermann et al. 

2011) has also recently published new results on the diffuse 
emission from the Cygnus region. An extended source, the so- 
called Cocoon, overlapping with MGRO J2031+41, has been 
found and its emission has been explained by freshly accelerated 
cosmic rays, trapped in a shell of photon-dominated emission 
formed by stellar winds and supernovae, as shown by mid- 
infrared observations. The spectrum from this region is hard, 
extending up to 100 GeV with no evidence of softening, and 
could be explained as the result of pulsar-accelerated particles 
or as an active super-bubble. The average diffuse emission from 
the region has also been analyzed at MeV to GeV energies (The 
Fermi LAT Collaboration 2012a). Despite the very rich source 
population, this emission is similar to that of the local interstellar 
space, and there is no evidence that it is necessary to include an 
extra contribution in the model, resembling the diffuse excess 
previously measured by Milagro (Abdo et al. 2007a, 2008). Most 
recently, the ARGO-YBJ collaboration presented the results of 
a survey of the Cygnus region in the energy range of 600 GeV 
to 7 TeV. MGRO J2031+41 is observed with a significance of 
6.4a and a flux consistent with previous Milagro results, but 
there is no significant detection of MGRO J2019+37 (Bartoli 
et al. 2012). 

Here, we present a new analysis of the last three years of 
data collected with the Milagro experiment (2005-2008). An 
improved gamma-hadron separation and a newly developed 
technique are exploited to reconstruct the energy spectra of 
gamma rays from the sources in the Cygnus region. The best 
fits for the spectra of the two brightest Milagro sources are 
presented: MGRO J2019+37, a source observed in the complete 
Milagro data set collected over eight years of operation with a 
pre-trials significance in excess of 12a between 1 and 100 TeV, 
and MGRO J2031+41, a source observed in the complete data 
set with a pre-trials significance in excess of la (Abdo et al. 
2009b). We compare our spectra with results from the HEGRA, 
MAGIC, Whipple, and ARGO-YBJ experiments (Aharonian 
et al. 2005; Albert et al. 2008; Lang et al. 2004; Bartoli et al. 

2012) . 


2. ANALYSIS TECHNIQUE 

The Milagro detector was located in the Jemez Mountains in 
New Mexico, at an altitude of 2630 m a.s.l. It was operated from 
2001 to 2008, and, in its final configuration, consisted of two 
components: (1) a central pond (60 x 80 m 2 , 8 m deep) with two 
photomultiplier tube (PMT) layers, a shallow air shower layer 


consisting of 450 PMTs, and a deep muon layer consisting of 
273 PMTs and (2) an array of 175 single-PMT tank detectors 
surrounding the pond covering an area of 40,000 m 2 (Atkins 
et al. 2003). 

A detailed description of the analysis method and the 
parameters used here can be found in the Milagro paper on the 
spectral measurement of the Crab Nebula (Abdo et al. 2012). 
The method has been confirmed to be applicable for cosmic- and 
gamma-ray energies between 1 and 100 TeV. Using cosmic-ray 
data, systematic uncertainties of this method are estimated to be 
~30% on the overall flux and ~0. 1 on the spectral index (see 
Section 4 of Abdo et al. 2012 for more details). 

The new element in Abdo et al. (2012) and this paper, with 
respect to the previous approach (Abdo et al. 2007b, 2009b), is 
the introduction of an estimator, T, for the energy of the primary 
particle initiating the extensive air shower, based on the number 
of PMTs hit for each event. The T parameter is defined as the 
sum of two fractions: 


ZVa.s. At. a. 

iwJive yylive ’ 
"A.S. iv T.A. 


( 1 ) 


where ZVa.s. and ZV S are the number of PMTs detecting the event, 
while and V]'(( are the number of functional PMTs in the 
air shower layer and in the tank array, respectively. Because 
the typical energy resolution of Milagro in each T bin is quite 
broad (50%-100%), the energy distributions between T bins 
overlap, and the median energy associated with a given T value 
is dependent on the spectral assumption for the observed source, 
the fit is performed in the measured T space rather than in the 
energy space. 

The reconstructed Milagro data contain information about 
the direction and T value of air shower events. According to 
their T value, events are sorted into nine bins (0.2 ^ T ^ 2), 
resulting in a total of nine signal skymaps, binned in Oi 1 x 0? 1 
pixels. For each T bin, background maps are calculated using 
the Direct Integration method with the typical 2 hr integration 
duration (Atkins et al. 2003; Abdo et al. 2012). The two 
brightest regions of interest in the sky, a 2° x 2° area around 
the position of the Crab Nebula and a band around the Galactic 
plane (— 2?5 < b < 2:5), are excluded when calculating the 
background (Abdo et al. 2012). Rather than discriminating 
between gamma-ray and cosmic-ray initiated air showers, a 
weight is applied to all measured events, in both signal and 
background maps, where gamma-ray-like events are assigned a 
higher weight than cosmic-ray-like events (Abdo et al. 2012). 
The gamma-ray excess with respect to the estimated background 
is calculated as the difference between signal weights and 
background weights, and is filled into the so-called excess map. 
Based on these excess maps, we compute the T distributions 
used in the energy fit. 

Excess maps are smoothed according to the detector angular 
resolution (or point-spread function, PSF), which is a function of 
T and varies between 0: 3 and 0. 7. The measured T distribution 
for the target source in a 0?1 x Oil bin, as a result, is the 
average excess coming from a PSF-wide region around the 
nominal source position. Since the background is measured 
with the direct integration method removing the Galactic plane 
( b > 2)5 and b < —2)5), the measured excess from sources 
in the Galactic plane includes any other diffuse or extended 
emission possibly present in the vicinity of the source. The 
Galactic diffuse background was estimated to contribute up to 
15% of the flux at 35 TeV for the weakest Galactic BSL source 
(Abdo et al. 2009b). 
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Table 1 

Source Positions 
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Survey 

R.A. 

(deg) 

Decl. 

(deg) 

1 o Error Ellipse (Semi) 
(deg) 

Width (cr) 
(deg) 

MGRO J2019+37 

Two-dimensional Gaussian 3 

304.63 

36.88 

0.11,0.08 

0.70 

Fermi position (2FGL J2021. 0+365 l) c 

305.2702 

36.8634 

0.0081,0.0078 


Milagro 2007 d 

304.98 

36.66 

0.20,0.13 

l.l d 

MGRO J203 1+41 

Two-dimensional Gaussian 3,b 

307.18 

41.31 

0.42, 0.26 

1.8 

Fermi position (2FGL J2032.2+4126) c 

308.0622 

41.4369 

0.0142,0.0133 


Milagro 2007 d 

307.98 

41.51 

0.52, 0.33 

3.0 d 


Notes. Positions of MGRO J2019+37 (top) and MGRO J203 1+41 (bottom) for the two-dimensional Gaussian fit, Fermi position of 
associated sources and previous Milagro survey. 

a The table lists statistical errors only. The systematic pointing error is <0?3. 

b Studies with subsets of data result in an additional systematic uncertainty of <0?5 in the case of MGRO J2031+41. 
c Abdo et al. (2009b). 

d The number reported is the diameter of a top-hat function that best fits the angular extent of the source (Abdo et al. 2007b). 


We test the emission from the two sources for two different 
spectral hypotheses: a power law, defined as 

d.N ( E \““ 

~dE ~ Nl0TeV V 10 TeV ) ’ 

where /VkitcV is the normalization scale set at 10 TeV and a is 
the spectral index, and a power law with an exponential cutoff, 
defined as 


dN / E \~ a / E\ 

— - N WTeV (^ 1()XeV J ex P(-^J’ () 

where E c is the cutoff energy. 

The fit to the data is performed comparing the measured ex- 
cess to simulations. First, a set of simulated data is generated 
varying sensitive parameters (i.e., A/foTeV, a, and E c ). The best 
spectral parameters and the corresponding fit probability are 
then found using a / 2 minimization, comparing the measured 
and the simulated T distributions. Uncertainties on fit param- 
eters are computed using la contours of x 2 histograms (as 
discussed in Sections 3.1 and 3.2). The uncertainty is defined as 
the distance between the best-fit value and the lower and upper 
edges of the la contour. 
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Figure 1. MGRO J20 19+37: distributions of the parameter used to estimate 
the photon spectrum, T . The spectrum fit is performed in the T space. Black 
squares represent the data. Red triangles and blue circles represent the simulated 
distributions for the best fit assuming a power law with a cutoff and a simple 
power law, respectively. The unit on the y-axis is the weighted excess per day, the 
unit on the x-axis is the T value. The error bars are the statistical uncertainties. 
(A color version of this figure is available in the online journal.) 


3. RESULTS AND DISCUSSION 

The results presented here are obtained from the analysis of 
the Milagro data taken with the detector in its final configuration. 
The start date is 2005 October 22 20:39:16 GMT and the stop 
date is 2008 April 15 00:02:53 GMT, corresponding to a total of 
906 days of observation, resulting in 832 integrated days after 
data quality cuts. 

The positions of MGRO J20 19+37 and MGRO J2031+41 are 
obtained fitting Milagro excess maps with a two-dimensional 
Gaussian function with equal widths in right ascension and 
declination (where the right ascension is corrected to get the 
true spherical distance). The fit is performed in a square region 
of ±3° around the excess. Table 1 shows the positions used 
in this paper and a comparison to other surveys. Table 2 lists 
the total gamma-ray event excess in each T bin for the two 
sources. 


3.1. MGRO J2019+37 Spectrum 

In Figure 1, the T distribution of data is shown with the 
two simulated T distributions for the best-fitting power law and 
power law with exponential cutoff models. These distributions 
show the weighted excess coming from the source position after 
background subtraction. The excess from MGRO J20 19+37 is 
compared to the excess measured from the entire Cygnus region 
(65° < l < 85° and —2° < b < 2°). We find that the emission 
from MGRO J20 19+37 accounts for ~12% of the excess from 
the entire region. 

Figures 2 and 3 show the la and 2a contours from the x 2 fit 
for the two tested hypotheses, la contours are used to compute 
uncertainty on the parameters as described in Section 2. 

The best-fit spectra for the two tested hypotheses 
(Equations (2) and (3)) are shown in Figure 4. The best-fit 
spectral parameters are summarized in Table 3. 
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a 

Figure 2. MGRO J2019+37: ltr and 2 a CL contours for the power-law fit. The 
point indicates the best-fit result (minimum x 2 )- The x 2 increments for the lcr 
and 2 a contours are 2.30 and 6.18, respectively. 

(A color version of this figure is available in the online journal.) 


With the power-law hypothesis, the best-fit parameters for the 
normalization and the spectral index are /VioTeV = (4.1^jg) x 
10 -10 s -1 m -2 TeV -1 and a = 2.78^%, respectively. The y 2 
is 16.12 for 7 degrees of freedom (dof), which gives a y 2 fit 
probability of 0.024. 

In the case of the power law with cutoff hypothesis, the 
best-fit parameters are IVioTeV = 7^ 5 2 x 10 -10 s -1 m -2 TeV -1 , 
a = 2.tP°j 5 0 , and E c = 29+ 5 ° 6 TeV. The y 2 is 1.924 (6 dof), 
which gives a y 2 fit probability of 0.93. 

The y 1 fit probabilities suggest that the power law with 
a cutoff model (/ 2 probability = 93%) fits the data better 
than a simple power law (/ 2 probability = 2.4%). An /-’-test 
(Bevington & Robinson 2003; Bates & Watts 2007) for the 
drop in x 2 when adding the cutoff parameter to the fit gives a 
probability of 6 x 10 -4 of having seen a larger change, if the 
simpler power law without exponential cutoff were the correct 
model. To test the robustness of the result, since the highest 
T bins may be subject to larger systematic uncertainties, we 
repeat the spectral fit excluding the last, and the last two bins of 
the distribution. The extracted parameter values do not change 
significantly, and despite the inclusion of fewer data, the simple 
power law remains disfavored. The /-’-test for the drop in y 1 
when adding the cutoff parameter to the spectral fit in only 8 (7) 



Figure 4. MGRO J20 19+37: energy spectra. The best fit is obtained for a power 
law with cutoff model (in red). The power-law model is also shown (in blue). 
The shadowed area represents the la band, obtained by varying the parameters 
within the la contour. ARGO-YBJ 90% CL upper limits for MGRO J20 19+37 
are shown in black (Bartoli et al. 2012). 

(A color version of this figure is available in the online journal.) 

Thins gives a probability of 4 x 10 -3 (6 x 10 -2 ) of having seen 
a larger change, if the simpler power law without exponential 
cutoff were the correct model. 

Previous Milagro analyses quoted the flux of MGRO 
J20 19+37 at 20 and 35 TeV, respectively (Abdo et al. 2007b, 
2009b). Those values, using a different analysis technique, are 
in agreement with the results presented here. An independent 
analysis (Allen 2007), which used a different energy estimator 
for the particle initiating the air shower and a different param- 
eter to distinguish between gamma and cosmic rays, is also 
consistent with our results. ARGO-YBJ (Bartoli et al. 2012) has 
not observed any significant emission from MGRO J2019+37 
in the range 600 GeV to 7 TeV. The explanations put forward 
in (Bartoli et al. 2012) are an insufficient exposure above 5 
TeV or a possible time variability of MGRO J2019+37. Due to 
the extension of MGRO J2019+37, the latter explanation is not 
without problems. The 90% confidence level (CL) upper limit 
from ARGO-YBJ is consistent with our best-fitting model, a 
power law with an exponential cutoff (see Figure 4). On the 



Figure 3. MGRO J2019+37: la and 2 a CL contours for the power law with cutoff fit. Left panel: AjoTeV vs. a projection. Right panel: E c vs. a projection. The points 
indicate the best-fit result (minimum x 2 )- The x 2 increments for the la and 2a contours are 3.53 and 8.02, respectively. 

(A color version of this figure is available in the online journal.) 
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Table 2 

Weighted Gamma-ray Event Excess 


Abdo et al. 


T 

0.2-0.4 

0.4-0.6 

0.6-0.8 

0.8-1 .0 

1.0-1. 2 

1.2-1. 4 

1.4-1 .6 

1.6-1. 8 

1. 8-2.0 

MGRO J2019+37 

1 1 7 +57 
iA '-57 

079+90 

• 3/z -90 

39329! 

5862*104 

A zq+12 0 
468_i2o 

8362*137 

1 09 5 +208 
1 UVD— 208 

2962 9 9 °o 

_4+14 

^-14 

MGRO J203 1+41 

218^57 

216 +90 
z 10—90 

28421 

2772*106 

901+122 

zyi -122 

1 622*143 

9+214 
z — 214 

100 -86 

_4+ll 

^-11 


Notes. Total weighted gamma-ray event excess in each T bin in 832 integrated days of data. 


Table 3 
Fit Results 


Source 

Mo TeV (10 -10 s -1 m- 2 TeV -1 ) 

a 

E c (TeV) 

X 2 /dof 

Prob. 

MGRO J20 19+37 

4 1+0-6 

9 7 Q+ 0.10 
'°- 0.10 

OO 

16.12/7 

0.024 

MGRO J20 19+37 

7+5 

—2 

2 0 +0 - 5 
z - u -i.o 

2Q+50 

zy -16 

1.924/6 

0.93 

MGRO J203 1+41 

9 1 + 0.6 
Z ‘ — 0.6 

9 99 +O .23 
J - zz — 0.18 

OO 

5.059/7 

0.65 

MGRO J203 1+41 

4T + 157 

9 7 + 0.7 
Z- ' —3.3 

21+ 00 
Z — 18 

3.352/6 

0.76 


Notes. Best-fit results for MGRO J2019+37 and MGRO J2031+41 using both the power law and the power law with cutoff hypotheses. 
The table lists statistical errors only, the systematic flux error is ~30% and the systematic error on the spectral index is ~0.1. 



7 

Figure 5. MGRO J2031+41: distributions of the parameter used to estimate 
the photon spectrum, T . The spectrum fit is performed in the F space. Black 
squares represent the data. Red triangles and blue circles represent the simulated 
distributions for the best fit assuming a power law with a cutoff and a simple 
power law, respectively. The unit on the y-axis is the weighted excess per day, the 
unit on the x-axis is the T value. The error bars are the statistical uncertainties. 
(A color version of this figure is available in the online journal.) 

other hand, the simple power-law model (a = 2.782o\o), al- 
ready disfavored by the E-test, does not agree with the ARGO- 
YBJ results. As mentioned in the Introduction, the VERITAS 
count map of the region shows an excess associated with MGRO 
J2019+37 and a separate source in the vicinity, VER J2016+372, 
for which VERITAS reported a flux of 1% of the Crab Nebula 
flux above 1 TeV and a hard spectrum with a photon index 
of 2.1+° 0 4 4 (Aliu et al. 2011). The hard spectrum VERITAS re- 
ports for VER J20 16+372 agrees well with the hard spectrum 
measured by Milagro for MGRO J2019+37 in the case of the 
power law plus exponential cutoff assumption. However, if VER 
J201 6+372 contributes to the flux of MGRO J2019+37, then it is 
a very small amount, since the flux of MGRO J2019+37 is about 
80% of the Crab Nebula and VER J2016+372 is separated from 
the larger VERITAS excess associated with MGRO 12019+37. 
VERITAS has not reported a spectrum and flux for the latter 
excess, yet. 



a 


Figure 6. MGRO J203 1+41 : la and 2a CL contours for the power-law tit. The 
point indicates the best-fit result (minimum y 2 ). 

(A color version of this figure is available in the online journal.) 


3.2. MGRO J203 1+41 Spectrum 

The T distribution from data and those simulated for the 
two best-fit models are shown in Figure 5. MGRO J2031+41 
accounts for ~6% of the excess from the entire Cygnus region. 

Figures 6 and 7 show the 1 a and 2er contours from the y_ 2 fit 
for the two tested hypotheses. 

Best-fit spectral parameters are summarized in Table 3. 

For the power-law hypothesis, the best-fit parameters are 
AlOTev = 2.H° 0 6 6 x 10 -1 ° s -1 m -2 TeV -1 , and a = 3.22+° 0 23 s . 
The x 2 is 5.06 (7 dof), with a x 2 fit probability of 0.65. 

For the power law with cutoff hypothesis, the best-fit param- 
eters are /VioTeV = 52.3 57 x 10 -10 s -1 m -2 TeV -1 , a — 2.7*33, 
and E c — 21*“ 8 TeV. In this case, we are not able to constrain 
the upper limit of the 1 a uncertainty for the cutoff energy (see 
the right panel of Figure 7). The maximum E c simulated value 
is 1000 TeV, an order of magnitude above the highest detected 
energy. The x 2 is 3.35 (6 dof), with a x 2 fit probability of 0.76. 
The best-fit spectra for the two models are shown in Figure 8. 

The improvement in the fit obtained with the more complex 
model is not significant, according to the E-test, and none of the 
parameters in the cutoff model are well constrained. Therefore, 
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Figure 7. MGRO J2031+41 : l<x and 2 a CL contours for the power law with cutoff fit. Left panel: AhoTeV vs. a projection. Right panel: E c vs. a projection. The points 
indicate the best-fit result (minimum x 2 )- The right panel shows that the present analysis is unable to constrain the upper limit of the energy cutoff within la. 

(A color version of this figure is available in the online journal.) 


Power Law 

Power Law with Cutoff 



Figure 8. MGRO J2031+41: energy spectra. The power-law model is shown in 
blue and the power law with cutoff model is shown in red. These two hypotheses 
give a similar x 2 fit probability. The pink cross is the Whipple flux at 0.6 TeV 
(Lang et al. 2004). The fine-dashed cyan line and the dashed yellow line are 
the HEGRA and MAGIC best fits (Aharonian et al. 2005; Albert et al. 2008), 
respectively. The dot-dashed black line is the ARGO-YBJ best fit (Bartoli et al. 
2012). The shadowed area represents the lcr band. 

(A color version of this figure is available in the online journal.) 

the simple power law is to be preferred over the power law with 
the exponential cutoff. 

The flux of MGRO J2031+41 at 20 and 35 TeV, measured 
by previous Milagro analyses (Abdo et al. 2007b, 2009b), is 
in agreement with the results presented here. Results from 
ARGO-YBJ, MAGIC, HEGRA, and Whipple (Bartoli et al. 
2012; Albert et al. 2008; Aharonian et al. 2005; Lang et al. 
2004) are shown in Figure 8. The HEGRA and MAGIC mod- 
els are mutually consistent, but they disagree with our best fit 
in terms of integral flux in the overlapping energy range. Be- 
tween 1 and 10 TeV, the flux measured by MAGIC accounts for 
only ~3% of the flux measured by Milagro. The spectrum as 
measured by the air Cherenkov telescopes (ACTs) is also much 
harder (a ~2) than the Milagro power-law best fit (a ~3). 
This discrepancy could be explained by the following two facts. 


First, the angular resolution of HEGRA and MAGIC ( <0: 1 ) is 
much better than that of Milagro (0. 3-0:7). The same is true for 
ARGO-YBJ, which, with an angular resolution between 0.47 
and 2? 8, measured an emission consistent with the results pre- 
sented here. The Whipple flux was measured at 0.6 TeV, with 
a PSF of 0: 2 1 , and lies between the ARGO-YBJ measurement 
and the HEGRA and MAGIC measurements. The measurement 
presented in this paper would therefore include photons coming 
from a larger region around the nominal source position com- 
pared to ACTs. Second, the way the background is computed is 
also different. Whipple used the ON/OFF method and observed 
the source with a significance of 3.3cr, while MAGIC primarily 
was operated in wobble mode (5.6er). A possible explanation of 
the presented results is that TeV J2032+4 1 30, whose extension is 
slightly larger than the HEGRA and MAGIC angular resolution, 
is surrounded by an extended emission, ~ 4° x 3° according to 
the Milagro map (see Figure 9 discussed in Section 3.3), and ac- 
cording to the previously reported diameter of MGRO J203 1 +4 1 
of 3°0 that was derived using a two-dimensional Gaussian fit 
(Abdo et al. 2007b; see also Table 1). Therefore, we believe 
that Milagro and ARGO-YBJ are not able to disentangle the 
extended emission from the central source and observe a higher 
flux. MAGIC and HEGRA, on the other hand, using ON/OFF 
and wobble mode might count the extended emission as back- 
ground. As a result, the ACTs measure a fainter emission and a 
harder spectrum for TeV J2032+4130. 

3.3. Morphology 

The significance map of the Cygnus region for the last three 
years of Milagro data (2005-2008) is shown in Figure 9. An 
energy-dependent PSF smoothing is applied. The top panel 
shows the entire region in Galactic coordinates. Milagro po- 
sitions for MGRO J20 19+37 and MGRO J2031+41 are shown, 
as well as the high-significance Fermi LAT sources (> 10cr) from 
the BSL catalog as reported in the second Fermi catalog (Nolan 
et al. 2012) and the extended Fermi Cocoon (Ackermann et al. 
2011 ). 

The Fermi counts map for photon energies above 10 GeV 
overlaid with the Milagro significance contours in the Cygnus 
region is shown in Figure 10. Milagro detects an excess 
consistent with the position of the two pulsars PSR J2021+365 1 
and PSR J2032+4127, associated with MGRO J2019+37 and 
J2031+41, respectively, but no detection of significant emission 
coincident with PSR J202 1+4026. 
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Figure 9. PSF-smoothed significance maps of the Cygnus region. Top panel: black triangles mark the positions of the two Milagro sources with the la radial width 
(solid circles), white-filled circles mark the Fermi BSL source positions as reported in the second Fermi catalog (the lcr error is smaller than the marker size). The Fermi 
Cocoon position is marked by the cross and the dashed circle is its la width. Bottom left panel: zoomed MGRO J20 19+37 significance map. The circle represents the 
la width. The positions of VER J2016+372 and 2FGL J2021. 0+3651 are also shown (see the text for a more detailed discussion of this region). Bottom right panel: 
zoomed MGRO J2031+41 significance map. The Milagro position with the la width is shown in black. The crosses mark the Whipple (green), HEGRA (magenta), 
and MAGIC (blue) positions with their uncertainty. 

(A color version of this figure is available in the online journal.) 



Figure 10. Fermi LAT counts maps of the Cygnus region above 10 GeV. We used Pass 7 photon data from 2008 August 4 15:43:37 to 2011 December 6 10:10:13 
GMT LAT Data Server, zenith angle <107°, data_qual=1, lat_config= 1 . A Gaussian smoothing is applied for display purposes. Circles mark the position of the 
Fermi BSL sources as reported in the second Fermi catalog. The green, blue, and black lines are the Milagro 3a, 5a, and 1 la level significance contours, respectively. 
(A color version of this figure is available in the online journal.) 


MGRO J2019+37 is observed with a pre-trials significance 
of 12.4cr in the complete Milagro data set collected over eight 
years of operation. Its extension is obtained with the two- 
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dimensional Gaussian fit discussed in Section 3. The fit results 
in a = 0:7, slightly larger than the angular resolution of 
the detector in its final configuration. There is evidence from 
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VERITAS that the extended source Milagro observes is the 
result of the superimposition of close point-like sources (Aliu 
et al. 2011). 

MGRO J2031+41, observed with a pre-trials significance of 
7.6er in the complete Milagro data set, is an extended source, 
with a = 178. It shows a central core with a higher significance 
(>6cr pre-trials) surrounded by a broader region (approximately 
3° x 4°) with a lower significance (>4a pre-trials). TeV 
J2032+4130 lies within the central core, and the Fermi Cocoon 
la circle (a = 2°) completely includes MGRO J2031+41. 
However, the Fermi analysis (Ackermann et al. 2011) shows 
that the Cocoon has an elongated shape, poorly correlating with 
the Milagro contours, and the possible association between the 
Cocoon and the extended Milagro excess remains unclear. 

The spectral fit results do not change if the maximum 
significance position is chosen instead of the two-dimensional 
Gaussian fit position. 

4. CONCLUSIONS 

We present the spectra of the two brightest Milagro sources 
in the Cygnus region using a new analysis technique applied to 
the last three years of data collected by the Milagro experiment. 

MGRO J2019+37 is observed with a pre-trials significance 
over 12(7 in the complete Milagro data set during eight years 
of operation. Its emission is well fitted by a power law with 
an exponential cutoff (E c = 29* 3 g TeV) and a hard asymp- 
totic spectral index (a = 2.0 +(i [ 'q). The simple power law is 
disfavored. An /-’-test (Bevington & Robinson 2003; Bates & 
Watts 2007) for the drop in x 2 when adding the cutoff param- 
eter to the fit gives a probability of 6 x 10 4 of having seen 
a larger change, if the simpler power law without exponential 
cutoff were the correct model. When excluding the highest (two 
highest) T bin(s), this value goes up to 4 x 10~ 3 (6 x 10 -2 ). 
The TeV excess measured by Milagro from MGRO J2019+37, 
spatially associated with the Fermi LAT pulsars J20 18.0+3626 
and J202 1.0+3651, has been confirmed by VERITAS, and it is 
likely produced by several nearby unresolved sources. The flux 
of these sources has not been reported by VERITAS, yet, though 
a separate source in the vicinity, VER J2016+372, is measured 
to have 1% of the Crab Nebula above 1 TeV and a spectrum 
with a photon index of 2 .Dq 4 (Aliu et al. 2011). ARGO-YBJ 
does not detect a significant emission from MGRO J2019+37, 
but the 90% CL upper limits do not conflict with the Milagro 
best-fitting model. 

The emission from MGRO J2031+41 (7.6a pre-trials signif- 
icance in the complete Milagro data set) is modeled by a power 
law with a = 3.22^o 23 8 . Our result, in particular the integral 
flux in the overlapping energy range, is consistent with previous 
measurements by ARGO-YBJ, but it disagrees with Whipple, 


HEGRA, and MAGIC results for TeV J2032+4130, most likely 
because of the different PSF of the instruments and different 
background subtraction methods. MGRO J2031+41 appears to 
show a structured morphology, produced by the superimposition 
of a central point-like source, coincident with TeV J2032+4130, 
and an extended emission that contributes to the overall spec- 
trum Milagro measures for MGRO J2031+41. This extended 
emission is possibly produced by either unresolved sources or 
interactions of cosmic rays with the local interstellar medium. 
The correlation between the TeV extended emission and the 
overlapping Fermi Cocoon is unclear and needs further studies. 

HAWC, the next generation water Cherenkov observatory, 
will be able to produce a more accurate analysis of the TeV 
emission from the Cygnus region, with its improved sensitivity 
(10-15 times better than Milagro) and angular resolution. 
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